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Abstract 

Evaluation of default correlation is an important task in credit risk analysis. In 
many practical situations, it concerns the joint defaults of several correlated firms, 
the task that is reducible to a first passage time (FPT) problem. This task represents 
a great challenge for jump-diffusion processes (JDP), where except for very basic 
cases, there are no analytical solutions for such problems. In this contribution, we 
generalize our previous fast Monte-Carlo method (non-correlated jump-diffusion 
cases) for multivariate (and correlated) jump-diffusion processes. This general- 
ization allows us, among other things, to evaluate the default events of several 
correlated assets based on a set of empirical data. The developed technique is an 
efficient tool for a number of other applications, including credit risk and option 
pricing. 
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1 Introduction 

In the financial world, individual companies are usually linked together via economic 
conditions, so default correlation, defined as the risk of multiple companies' default 
together, has been an important area of research in credit analysis with applications to 
joint default, credit derivatives, asset pricing and risk management. 

Nevertheless, the development of efficient computational tools for modeling default 
correlation is lagging behind its practical needs. Currently, there are two dominant 
groups of theoretical models used in default correlation. One is a reduced form model, 
such as in ifTZl that uses a Copula function to parameterize the default correlation. Re- 
cently, Chen et al J4] have translated the joint default probability into bivariate normal 
probability function. 
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The second group of models for default correlation is a structural form model. Zhou 
lfl6l and Hull et al [8 1 are the first to incorporate default correlation into the Black-Cox 
first passage structural model. They got similar closed form solutions for two assets. 
However, their models cannot easily be extended to more than two assets. Furthermore, 
the developed models do not include jump-diffusion processes. 

As demonstrated in flTI, jump risk becomes an important factor in credit risk anal- 
ysis. It is now widely acknowledged that the standard Brownian motion model for 
market behavior falls short of explaining empirical observations of market returns and 
their underlying derivative prices ifTTI . The (multivariate) jump-diffusion model that 
provides a convenient framework for investigating default correlation with jumps be- 
comes more readily accepted in the financial world. 

One of the major problems in default analysis is to determine when a default will 
occur within a given time horizon, or in other words, what the default rate is during 
such a time horizon. This problem is reduced to a first passage time (FPT) problem 
that can be formalized on a basis of a certain stochastic differential equations (SDE). 
It concerns the estimation of the probability density of the time for a random process 
to cross a specified threshold level. Unfortunately, after including jumps, only special 
cases have analytical solutions. For most practical cases, closed form solutions are 
unavailable and we can only turn to the numerical procedures. 

Monte Carlo simulation is such a candidate for solving the SDE, arising in the 
context of FPT problem. In conventional Monte Carlo methods, we need to discretize 
the time horizon into small enough intervals in order to avoid discretization bias iflOl . 
and we need to evaluate the processes at each discretized time which is very time- 
consuming. Many researchers have contributed to the field and enhanced the effi- 
ciency of Monte Carlo simulation. Atiya and Metwally 0] [T3ll have developed a fast 
Monte Carlo-type numerical method to solve the FPT problem. Recently, we have re- 
ported an extension of this fast Monte Carlo-type method in the context of multiple 
non-correlated jump-diffusion processes Q3). 

In this contribution, we develop a methodology for solving the FPT problem in the 
context of multivariate jump-diffusion processes. In particular, we have generalized the 
reported fast Monte-Carlo method (non-correlated jump-diffusion cases) for multivari- 
ate (and correlated) jump-diffusion processes. The paper is organized in the following 
way: Section [2] details the description of our model. The algorithms are described in 
Section [3] Section [4] contains the simulation results and discussions, followed by the 
conclusions given in Section[5] 

2 Developing multi-dimensional model 

2.1 Default correlation 

In the market economy, individual companies are inevitably linked together via dynam- 
ically changing economic conditions. Take two firms A and B as an example, whose 
probabilities of default are P4 and Pg, respectively. Then the default correlation can 
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be defined as 



PAB 



Pab - PaPb 



(1) 



^Pa(1-Pa)Pb(1-PbY 
where Pab is the probability of joint default. 

Now, we can write Pab as Pab = PaPb + pab\JPa(Jl - Pa)Pb{^ - Pb)- If 
we assume Pa = Pb = p « 1, then we have Pab ~ P 2 + PabP ~ PabP- Thus 
it is apparent that the default correlation pab plays a key role in the joint default with 
important implication in the field of credit analysis. Zhou |[T6l and Hull et al J8) were 
the first to incorporate default correlation into the Black-Cox first passage structural 
model. 

In |[T6l Zhou has proposed a first passage time model to describe default correla- 
tions of two firms under the "bivariate diffusion process": 



(2) 



where pi and p 2 are constant drift terms, Z\ and z 2 are two independent standard 
Brownian motions, and O is a constant 2x2 matrix such that 
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The coefficient p reflects the correlation between the movements in the asset values 
of the two firms. Based on Eq. (O, Zhou has deduced the closed form solution of 
default correlations of two assets lfl6l . However, none of the above models has included 
possible jumps. Apparently, jumps have more significant importance in the default 
correlation than often perceived. Indeed, simultaneous jumps may enhance the chance 
of simultaneous defaults which increases the correlation defaults. 



2.2 Multivariate jump-diffusion process 

Let us consider a more general case. In a complete probability space (Q, F, P) with 
information filtration (-Ft). Suppose that X t = ln(Vt) is a Markov process in some 
state space D C W l , solving the stochastic differential equation (7] 

dX t = p(X t )dt + a{X t )dW t + dZ t , (3) 

where W is an (F t )-standard Brownian motion in R"; p, : D — > M", a : D — > M nxra , 
and Z is a pure jump process whose jumps have a fixed probability distribution v on 
K™ such that they arrive with intensity {X(X t ) : t > 0}, for some A : D — > [0, oo). 
Under these conditions, the above model is reduced to an affine model if Q: 

p{X u t) =K + K 1 X t 

(a(X t ,t)a(X t ,t) T ) ij = (Ho)a + (H^X, 

\(X t ) =h + h- X t , (4) 

where K = (K Q ,Ki} G K™ x K" x ™, H = (Hq,^) g M" x " x R™x»x«, I = 
(l ,h) G R" x R" x ™. 
If we assume that, 
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1 . Each Wt in Eq. ([3]l is independent; 

2. K\ — 0, Hi — and l\ — in (4) that means the drift term, the diffusion 
process (Brownian motion) and the arrival intensity are independent with state 
vector X t ; 

3. The distribution of jump-size Z t is also independent with respect to X t . 
In this scenario, we can rewrite Eq. as 

dX t = iidt + adW t + dZ t , (5) 

where 

/i = Ko, (7(j T = Ho, X = Iq. 

At first sight, Eq. <(5j is similar to Eq. ©, but Eq. (O is a much more general model 
that can be applied to multiple firms, and where jumps have been taken into account. 

2.3 First passage time distribution 

Let us consider a firm i, as described by Eq. (|5j, such that its state vector Xi satisfies 
the following SDE: 

dXi = mdt + *^2 CijdWj + dZi = ^dt + a t dWi + dZ. h (6) 

3 

where Wi is a standard Brownian motion and Oi is: 




We assume that in the interval [0, T], total number of jumps for firm i is Mi times 
of jumps. Let the jump instants be Ti, T2, ■ ■ ■ , Tm^ Let To = and = T. The 

quantities tj equal to interjump times, which is Tj — Tj_\. Let Xi(T~) be the process 
value immediately before the jth jump, and Xi{T^~) be the process value immediately 
after the jth jump. The jump-size is Xi{T^) — Xi(T^), we can use such jump-sizes 
to generate Xi(T^) sequentially. 

In a structural model, a firm defaults when the firm assets value Vt falls below a 
threshold level Dv{t). In this contribution, we use an exponential form, defining the 
threshold level by Dy{t) = k exp(~ft) as proposed in 0161 , where 7 can be interpreted 
as the growth rate of firm's liabilities. Coefficient «, in front, captures the liability 
structure of the firm, which is usually defined as a firm's short-term liability plus 50% 
of the firm's long-term liability. As mentioned before, X t = ln(Vt), then the threshold 
of X t is D(t) = 7 < + 1ii(k). 

Atiya et al (TJ have deduced a one-dimensional first passage time distribution in 
time horizon [0, T]. In order to evaluate multi-firms, we obtain multi-dimensional for- 
mulas and reduce them to computable forms. 
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First, let us define A4 (t) as the event that process crossed the threshold level Di (t) 
for the first time in the interval [t,t + dt]. Then we have 

9ij (t) = P (A(t) G dtlX^T^), Xi{T7)). (7) 



If we only consider one interval [Tj_i, Tj], we can obtain 

XATf ,) - DAt) 
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/ [X^Tp - D, t (t) - ^(T, - t)} 

* cxp ( wrwf — ) 

* 6XP ( J' W 

where 

[X 4 (T+ 1 )-X J (Tr) +Al4Tj f\ 



2/j = 75= CX P 1 o 2 1 • 

After getting these results in one interval, we combine the results to obtain the 
density for the whole interval [0,T]. Let B(s) be a Brownian bridge in the interval 
[Tj-uTj] with BiTp) = XipTf^) and B(Tf) = X^Tf). Then the probability 
that the minimum of B(si) is always above the boundary level is 

Pii = p(„ B( Si )> D i (t)\B(T+_ 1 )=X i (Tf_ 1 ),B(T7)=X i (T j 




if JQ(Tr) > ^ (t) 



2[Xi (T+ x ) -£>; (*)] (T- ) (t)] 
^5 



otherwise. 



This implies that B(s\) is below the threshold level, which means the default hap- 
pens or already happened, and its probability is 1 — Py. Let L(si) = Li denote the 
index of the interjump period in which the time Sj falls in [T^-i, T^]. Also, let Ii 
represent the index of the first jump, which happened in the simulated jump instant: 

h = min(j : Xi(Tp > Di(t);k = l,...,j, and 

Xi{T+) > D t (t);k = 1 , and X t {T+) < D,(t)).(10) 

If no such Ii exists, then we set Ii = 0. 

By combining Eq. (HJ, © and ( fTol i. we get the probability of Xi crossing the 
boundary level in the whole interval [0, T] is 

P(A t ( Sl ) G ds|X ? ;(T+ J.XiCrr), j = 1, . . .,Mi + 1) 

giL, (si) rifeii 1 -fifc if ^» < ^ ° r 1% = °> 

5iL, rifcii 1 ^ + rifcii ^^(^ - t u ) if Lj = h, di) 

if > h, 

where 5 is the Dirac's delta function. 
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2.4 The kernel estimator 



For each firm, after generating a series of first passage times Sj, we use a kernel density 
estimator with Gaussian kernel to estimate the first passage time density (FPTD) /. As 
described in JT|, the kernel density estimator is based on centering a kernel function of 
a bandwidth as follows: 

- 1 N 

/=-^#(M-*i), (12) 

i=l 

where 

K(h,t-s i ) = - f ^expf (< ~ ^ 



The optimal bandwidth in the kernel function K can be calculated as fl4l : 

(/•oo \ —0.2 

2NV^ J (ft) 2 *) , (13) 

where TV is the number of generated points and f t is the true density. Here we use the 
approximation for the distribution, as a gamma distribution as proposed in (TJ: 



ft = T ^- ) t fj - 1 cxp(-at). (14) 
So the integral in Eq. ( fT3l becomes, 

r [r?dt f MSM (15) 

where 

W^i=A 2 , W 2 = 2AB, W3 = -B 2 + 2^4(7, W4 = 2_BC, W 5 = C 2 , 

and 

yl = a 2 , B=-2a(/?-l), C = (/3 - l)(/3 - 2). 

From Eq. ( fT31 l, follows that in order to get a nonzero bandwidth, we have constraint 
j3 to be at least equal to 3. 

After obtaining the estimated first passage time density /, the cumulative default 
rates can be written as, 

Pi(t) = I fi(r)dr. (16) 



3 Methodology of the solution 

In Section [2] we have reduced the solution of the original problem to a multivariate 
jump-diffusion model as described in Eq. ©. The first passage time distribution was 
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obtained in Section l231 As we have already noted, once jumps are included in the pro- 
cess, only for very basic applications closed form solutions are available lfTTl l2l. which 
in most practically interesting cases we have to resort to the numerical procedures. 

Let us recall the conventional Monte-Carlo procedure in application to the analysis 
of the evaluation of firm Xi within the time horizon [0, T], We divide the time horizon 
into n small intervals [0, fi], [£i,£2]> ■ ■ ■> [t n -i,T] as displayed in Fig. [TJ a )- In eacn 
Monte Carlo run, we need to calculate the value of X{ at each discretized time t. As 
usual in order to exclude discretization bias, the number n must be large. This pro- 
cedure exhibits substantial computational difficulties when applied to jump-diffusion 
process. 

Indeed, for a typical jump-diffusion process, as shown in Fig. [TJa)- let Tj-i and 
Tj be any successive jump instants, as described above. Then, in the conventional 
Monte Carlo method, although there is no jump occurring in the interval [Tj_i,Tj], 
yet we need to evaluate Xi at each discretized time t in [Tj-i,Tj]. This very time- 
consuming procedure results in serious shortcoming of the conventional Monte Carlo 
methodology. 



Conventional Monte-Carlo method 
J I I I 



h t 2 T M T } t n l T 

(a) 

Uniform sampling method 
j s, (uniform distribution) t 

J ~ l \ 1 1 J 

* ' 7 ' Brownian bridge » X ^ T 

(b) 



Figure 1 : Schematic diagram of (a) conventional Monte Carlo and (b) uniform sam- 
pling method. 

To remedy the situation, two modifications of the conventional procedure were re- 
cently proposed by Atiya and Metwally (T| Q~3l that allow us a potential speed-up of 
the conventional methodology in 10-30 times. One of the modifications, the uniform 
sampling method, involves sampling method. The other, inverse Gaussian density sam- 
pling, is based on the inverse Gaussian density method for sampling. Both methodolo- 
gies were developed for the univariate case. 
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In this article, we focus on the further development of uniform sampling (UNIF) 
method and extend it to multivariate jump-diffusion processes. The major improvement 
of the uniform sampling method is based on the fact that it only evaluates Xi at gener- 
ated jump instants while between each two jumps the process is a Brownian bridge (see 
Fig. |TJb)). Hence we just consider the probability of Xi defaults in (Tj-i, Tj) instead 
of evaluating Xi at each discretized time t. More precisely, in the uniform sampling 
method, we assume that the values of X^T+.i) and X t (Tp are known as two end 
points of the Brownian bridge, we generate a variable Sj with uniform distribution and 
by using Eq. (O we verify whether Xi(si) is smaller than the threshold level. If it is, 
then we have successfully generated a first passage time s, and can neglect the other 
intervals and perform another Monte Carlo run. 

In order to implement the UNIF method for our multivariate model as described in 
Eq. ([6), we need to consider several points as follows: 

1. Here, we focus on the firms rated in the same way, that is we assume that the 
arrival rate A for the Poisson jump process and the distribution of (Tj — 

are the same for each firm. As for the jump-size, we generate them by a given 
distribution, which can be different for different firms to reflect specifics of the 
jump process for each firm. 

2. We exemplify our description by considering an exponential distribution (mean 
value {It) for (Tj — Tj_i) and a normal distribution (mean value fij and stan- 
dard deviation erj) for the jump-size. We can use any other distribution when 
appropriate. 

3. An array IsDef ault (whose size is the number of firms denoted by A^ rm ) was 
used to indicate whether firm i has defaulted in this Monte Carlo run. If the firm 
defaults, then we set IsDefault(i) = 1, and will not evaluate it during this 
Monte Carlo run. 

4. Most importantly, in order to reflect the correlations of multiple firms, we need 
to generate correlated s;. As deduced in ll9l 1161 . the joint probability of firm 
i to default before Si and firm j to default before Sj satisfies the bivariate in- 
verse Gaussian distribution. We approximate the correlation between Si and Sj, 
denoted further by p Si . s as the default correlation of firm i and j. Then, we 
employ the sum-of-uniforms method [3 1 to generate correlated uniform numbers 

St- 

Our algorithm for multivariate jump-diffusion processes can be described as fol- 
lows. It is an extension of the one-dimensional case, described in |[Tl [T3l . 

Consider A^ rm firms in the interval [0, T]. First, we generate the jump instant Tj 
by generating interjump times (Tj — and set all the IsDefault(i) = to 

indicate that no firm defaults at first. 

From Fig. [TJb) and Eq. ©, we can conclude that, 

1. If no jump occurs, as described in Eq. ©, the interjump size (Xi(Tj~) — 
X^T^^) follows a normal distribution of mean Hi(Tj — Tj-i) and standard 
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deviation er, y/Tj — Tj~y. We get 

Tj^ + cnNiO^-Tj^) 

iV firm 

Tj-i)+ o-ifc^(0,r i -T i _ 1 ) > 
fe=l 

where the initial state is Ai(0) = Xj(Tg ). 

2. If jump occurs, we simulate the jump-size by a normal distribution or other dis- 
tribution when appropriate, and compute the postjump value: 

X i {Tf) = X i (T7) +Zi {T j ). 

After generating beforejump and postjump values Xj (Tj~ ) and Xi (Tj~ ) . As before, 
j = 1, M where M is the number of jumps for all the firms, we can compute Py 
according to Eq. To recur the first passage time density (FPTD) fi(t), we have to 
consider three possible cases that may occur for each non-default firm i: 

1. First passage happens inside the interval. We know if Xi(Tf_ t ) > Di(Tj-\) 

and Xi(T~) < Di(Tj), then the first passage happened in the time interval 
[Tj-i, Tj]. To evaluate when the first passage happened, we introduce a new 
viable fey, as fey = , so that + fey = Tj ~y^p^~^ ■ usm g the 

sum-of-uniforms method, we generate several correlated uniform numbers Sj in 
the interval of [Tj-i, Tj_i + fey], and if Sj also belongs to interval [Tj-i, Tj], 
then the first passage time occurred in this interval. We set IsDef ault(i) = 
1 to indicate firm i that has defaulted. Then, we can compute the conditional 
boundary crossing density gy(sj) according to Eq. (O. To get the density of 

the entire interval [0, T], we use = f ^Zp.T 1 J 9ij ( s i) * K(h opt , t — s,), 

where n is the iteration number of the Monte Carlo cycle. 

2. First passage does not happen in the current interval. If Sj doesn't belong 
to interval [Tj_i,Tj], then the first passage time has not yet occurred in this 
interval. 

3 . First passage happens at the right boundary of interval. If Xi (T+ ) < D^Tj) 

and Xi(Tj~) > Di{Tj) (see Eq. ([Toil), then T] t is the first passage time and 

Ii = j. We evaluate the density function using kernel function /i jTl (t) = 
K(h opt ,t- TjJ, and set IsDefault(z) = 1. 

Next, we increase j and examine the next interval and analyze the above three cases 
for each non-default firm again. After running times the Monte Carlo cycle, we get 
the FPTD of firm i as fi(t) = -k Y^ n =i fi,n{t)> as we U a s the cumulative default rates 

Pi(t)=Sofi(r)dr. 



Xi(Tf) ~ XiiT+J+iHpj- 



9 



4 Applications and discussion 



In this section, we will demonstrate how our model describes the default correlations 
of the firms, rated in the same way, via studying the historical data. We provide details 
on the calibration of the models we apply for this description. 



4.1 Default rates 

In Fig. [2 the black line of square is a set of historical default data of A-rated firrrQ 
taken from [16]. 

First, if we do not consider jumps, as assumed in Ifl6l , the firm defaults at time t 
with probability: 

m (-m^w) =2 .N (-*) , (17) 

where 



criVi J V y/i 



Xi(0)-ln(««) 
<Ji 

is the standardized distance of firm i to its default point and N(-) denotes the cumula- 
tive probability distribution function for a standard normal variable. 
If historical default rates are given, we can estimate Z± as follows: 

Zi = argmm ^ I — I , (18) 

where t) are the theoretical default probabilities (as determined by Eq. ( fT71 )) 

and Ai(i) are the historical default rates. For the A-rated firm considered here, the 
optimized Z^ value was evaluated in |[T6l as 8.06. By feeding the optimized Z^-value 
into Eq. ( fTTI i, we get the theoretical cumulative default rates without jumps, given in 
Fig. |2]by the line of circles. 

Now, let us consider the UNIF method, developed in Section [2~4l and [31 First, the 
developed Monte Carlo simulation allows us to obtain the estimated density /j (t) by 
using kernel estimator method. We get also the default rate Pj (t) for firm i. 

Then we minimize the difference between our model and the historical default data 
to obtain the optimized parameters in our model: 



argmin I 



Pi(tj)-Mtj) 



2^ 



4-^ \ t 



(19) 



j 



For convenience, we reduce the number of optimizing parameters by: 
1. Setting X(0) = 2 and 1h.(k) = 0. 



A-rated fiiTn stands for a specific kind of firm following the Moody's Investors Service's definition. 
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Figure 2: Historical, theoretical and simulated cumulative default rates. The theoret- 
ical value is calculated by using Eq. ( fTTI i. All the simulations were performed with 
Monte Carlo runs N = 100, 000, besides, for conventional Monte Carlo method, the 
discretization size of time horizon is A = 0.005. 
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2. Setting the growth rate 7 of debt value equivalent to the growth rate fx of the 
firm's value 1161 . so the default of firm is non-sensitive to /i. Setting /i = —0.001 
in our computation reported next. 

3. The interjump times (Tj — satisfy an exponential distribution with mean 
value equals to 1 . 

4. The arrival rate for jumps satisfies the Poisson distribution with intensity param- 
eter A, where the jump size is a normal distribution Z t ~ N{nz, crz)- 

As a result, we only need to optimize a, A, \Xz, &z f° r this firm, This is done by 
minimizing the differences between our simulated default rates and the given historical 
data. The minimization was performed by using quasi-Newton procedure implemented 
as a Scilab program. 

The optimized parameters for the A-rated firm are a = 0.09000984, A = 0.10001559, 
\iz = —0.20003641, and <rz = 0.50000485. Then, by using these optimized param- 
eters, we carried out a final simulation with Monte Carlo runs N = 100, 000. The 
simulated cumulative default rates by using the UNIF method are shown in Fig. [2] by 
the dash line. For comparison, we have carried out the conventional Monte Carlo sim- 
ulation with the same optimized parameters. The resulting simulated default rates are 
displayed by dotted line in Fig. [2] All the simulations reported here were carried out 
on a 2.4 GHz AMD Opteron(tm) Processor. The optimal bandwidth and CPU time are 
given in Table [T] 

Table 1: The optimal bandwidth h opt , and CPU time per Monte Carlo run of the sim- 
ulations. All the simulations were performed with Monte Carlo runs N = 100, 000; 
the discretization size of time horizon for the conventional Monte Carlo method was 
A = 0.005. 





Optimal bandwidth 


CPU time per Monte Carlo run 


Conventional Monte Carlo 


0.891077 


0.119668 


UNIF 


0.655522 


0.000621 



From Fig. [2] we can conclude that our simulations give similar results to the the- 
oretically predicted by Eq. (T% , and exceed them for short time horizon. The UNIF 
method gives exactly the same default curve as the conventional Monte Carlo method, 
but the former outperforms the latter substantially in terms of computational time. The 
UNIF methodology is much faster compared to the conventional method and is ex- 
tremely useful in practical application. 

4.2 Default correlations 

Our final example concerns with the default correlation of two A-rated firms (A,A). In 
Table |2] we provide the information on the default correlation of firms (A, A) for one-, 
two-, five- and ten-year. The values in the 2nd column were calculated by using the 
closed form solution derived in lfl6l . 

In order to implement the UNIF method, we use assumptions, similar to the ones 
before, in order to reduce the number of optimized parameters: 
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Table 2: Theoretical and simulated default correlations (%) of firms (A,A). The simu- 
lations were performed with Monte Carlo runs N = 100, 000. 



Year 


Ref. [16| 


UNIF 


1 


0.00 


0.00 


2 


0.02 


2.47 


5 


1.65 


6.58 


10 


7.75 


9.28 



1. Setting X(0) = 2 and ln(«;) = for all firms. 

2. Setting 7 = \i and /i = —0.001 for all firms. 

3. Since we are considering two same rated firms (A,A), we choose a as: 



where aa T = Hn such that 



aa =H = 



Oil 

021 



07 



CT12 

022 



P12C1CT2 



and 



a 



21 



a 



22 



P12 



PjU<721 + Q'120'22 
(7iCr2 



(20) 



(21) 



(22) 



In d22l . pi2 reflects the correlation of diffusion parts of the state vectors of the 
two firms. 

4. The arrival rate for jumps satisfies the Poisson distribution with intensity pa- 
rameter A for all firms, and we use the parameters optimized based on a single 
A-rated firm, i.e., A = 0.10001559 for all the firms. 

5. As before, we generate the same interjump times (Tj — that satisfies an 
exponential distribution with mean value equals 1 for all firms. Furthermore, the 
jump size is a normal distribution Zt ~ N(/j,Zi, ^z t ), and we use the parameters 
optimized from a single A-rated firm, i.e., [iz i = —0.20003641 and cr^ = 
0.50000485 for all the firms. 

As a result, there are only 4 parameters left to optimize: an, &12, 021 and 0-22- The 
optimization was carried out by using the quasi-Newton procedure implemented as a 
Scilab program. The resulting optimized parameters are an = 0.06963755, <Ti2 = 
0.02993134, a 2 \ = 0.03387809 and cr 2 2 = 0.06 69 1 001. We can easily get a 1 = 
0.0757976, a 2 = 0.0749978 and p 12 = 0.7673104. The parameter p u represents the 
correlation between diffusion parts of the state vectors of two firms. 

The simulated default correlations are displayed in the 3rd column of Table|2] Ob- 
serve that, the UNIF method gives a little larger default correlation compared to the 
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theoretical predicted by Eq. ( fTTb . This is mainly because our optimized px% is larger 
than 0.4 used in 1161 . and we have used the same interjump times (Tj — for all 

the firms. Nevertheless, the UNIF method gives the correct default correlation trend, 
as the default correlation becomes larger with increasing time. 

5 Conclusion 

We analyzed the first passage time problem in the context of multivariate and correlated 
jump-diffusion processes by extending the fast Monte Carlo-type numerical method - 
the UNIF method - to the multivariate case. We provided an application example of 
simulating default correlations confirming the validity of our model and the developed 
algorithm. Finally, we note that the developed methodology provides an efficient tool 
for further practical applications such as in credit analysis and barrier option pricing. 
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